#!/bin/bash -e

function info() {
echo Usage: `basename $0` [-i in.bed] [-r known_interval] in.bam
exit 1
}

while getopts  ":i:r:n:p:" opts
do
    case  $opts  in
        i) interval=$OPTARG;;
        r) realn_interval=$OPTARG;;
        p) out_prefix=$OPTARG;;
        n) data_thread_num=$OPTARG;;
        *) info;;
    esac
done
shift $(($OPTIND - 1))


if [ $# -lt 1 ]; then info; fi


. $var


vcf_path=${data_path}/vcf/gatk/

if test -z "$realn_interval"; then

if [ -n "$interval" ]; then

echo
echo
echo gatk RealignerTargetCreator recommand nt 24 mem 48
java $j_mem -jar $gatk \
	-T RealignerTargetCreator \
	-R $ref_genome \
	-I $1 \
	-o $out_prefix.realn.intervals \
	-log $out_prefix.realn.log \
	-nt $data_thread_num \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known $data_path/gatk/vcf/Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	-L $interval \
	-S SILENT -dt BY_SAMPLE -dcov 300 -l INFO
	
else

echo
echo
echo gatk RealignerTargetCreator recommand nt 24 mem 48
java $j_mem -jar $gatk \
	-T RealignerTargetCreator \
	-R $ref_genome \
	-I $1 \
	-o $out_prefix.realn.intervals \
	-log $out_prefix.realn.log \
	-nt $data_thread_num \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	-S SILENT -dt BY_SAMPLE -dcov 300 -l INFO
	
fi

realn_interval=$out_prefix.realn.intervals

else 
echo; echo; echo known interval: $realn_interval
fi


if test -z "$interval"; then

echo; echo; echo gatk IndelRealigner recommand sg 4 mem 4
java $j_mem -jar $gatk \
	-T IndelRealigner \
	-R $ref_genome \
	-targetIntervals $realn_interval \
	-I $1 \
	-o $out_prefix.realn.bam \
	-log $out_prefix.realn2.log \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	-S SILENT -dt BY_SAMPLE -dcov 300 -l INFO -LOD 0.4

else

# indel realign
echo; echo; echo gatk IndelRealigner recommand sg 4 mem 4
java $j_mem -jar $gatk \
	-T IndelRealigner \
	-R $ref_genome \
	-targetIntervals $realn_interval \
	-I $1 \
	-o $out_prefix.realn.bam \
	-log $out_prefix.realn2.log \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	-S SILENT -dt BY_SAMPLE -dcov 300 -l INFO -LOD 0.4 \
	-L $interval
	
fi	
	
	
. $cmd_done
exit 0
	
# LOD threshold above which the cleaner will clean
# This term is equivalent to "significance" - i.e. is the improvement significant enough to merit realignment? Note that this number should be adjusted based on your particular data set. For low coverage and/or when looking for indels with low allele frequency, this number should be smaller.
	
# -S --validation_strictness
# -dt --downsampling_type
# -l --logging_level / -l
# Set the minimum level of logging
# Setting INFO gets you INFO up to FATAL, setting ERROR gets you ERROR and FATAL level logging, and so on.
